GE ~ power log model from Ettema & Loras, 2009

library(tidyverse)
library(geomtextpath)
library(ggtext)
library(colorspace)
library(sysfonts)
library(showtext)

## use font & icons in ggplots
sysfonts::font_add_google(
    name = "Merriweather Sans", 
    family = "Merriweather Sans"
)
sysfonts::font_add(
    family = "fa-brands", 
    regular = r"(C:\Users\Jem\AppData\Local\Microsoft\Windows\Fonts\Font Awesome 7 Brands-Regular-400.otf)"
)
showtext::showtext_opts(dpi = 600)
showtext::showtext_auto()

## custom theme
theme_set(
    theme_bw(base_size = 12, base_family = "Merriweather Sans") +
        theme(
            plot.title = ggtext::element_textbox_simple(
                size = rel(1.2), lineheight = 1.1
            ),
            plot.subtitle = ggtext::element_textbox(
                colour = "grey35", face = "italic", hjust = 0, 
                lineheight = 1.1, margin = margin(t = 4, b = 4)
            ),
            plot.caption = ggtext::element_textbox(
                colour = "grey35", halign = 1
            ),
            panel.border = element_rect(linewidth = 0.8),
            axis.title = element_text(face = "bold"),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            legend.position = "none"
        )
)

social_caption <- "<span style='font-family:fa-brands'>&#xf09b; &#xe671;</span> jemarnold"

Post 1

Typical range for gross metabolic efficiency (GE) = 15-25%

GE logarithmically increases with power output, so higher fitness = higher power = higher GE

Most of the variance here comes from differences in cadence. The remaining variance comes from random cross-sectional individual differences

## digitized from Ettema & Loras, 2009. Efficiency in cycling: a review Figure 2.d
GE_data <- tribble(
    ~workload         , ~GE               ,
     24.1054592934376 ,  7.29508257346393 ,
     38.9202711580285 ,  8.62295161158726 ,
     31.8895188354463 , 10.4754096314559  ,
     79.5982326093705 , 12.9672142287399  ,
     48.2109185868752 , 12.9672142287399  ,
     48.96421695563   , 13.8688529856967  ,
     49.2153105114338 , 14.9508206497105  ,
     65.0345231512551 , 15.1147541286731  ,
     74.0740681736263 , 14.9508206497105  ,
    119.271810986825  , 14.7704928983192  ,
     94.9152581375841 , 15.5573775266028  ,
     94.1619597688293 , 15.9999997688668  ,
     63.7790376708926 , 16.8524591755342  ,
     97.1751355425052 , 17.836065827638   ,
     97.92843391126   , 17.5901644535284  ,
    100.439422573328  , 17.0983605496438  ,
     99.6861242045735 , 19.0983612430431  ,
    124.042694755158  , 15.7540995504229  ,
    129.315765635098  , 17.3442630794189  ,
    129.566841489559  , 17.9836073455031  ,
    127.809168897589  , 18.2131156028495  ,
    128.813543120804  , 18.5245905996771  ,
    125.047068978373  , 18.5245905996771  ,
    130.06904630251   , 19.0327887759906  ,
    144.883849316429  , 19.6885250031723  ,
    145.13497827492   , 19.2622958776714  ,
    144.632755760625  , 18.9180334916518  ,
    150.156920196369  , 19.3934431231077  ,
    149.654733084762  , 19.7049181199354  ,
    149.654733084762  , 19.9016401437555  ,
    150.156920196369  , 20.9344273018143  ,
    149.403639528958  , 21.7213119300979  ,
    149.905826640566  , 16.3606564273152  ,
    159.698687733035  , 15.7377052779942  ,
    180.288784141189  , 17.3606561961821  ,
    164.720647355828  , 18.3278697315227  ,
    159.949781288839  , 19.0000002311331  ,
    161.205284470545  , 19.0983612430431  ,
    163.967366688417  , 19.5081972517809  ,
    174.764584302758  , 18.8688529856967  ,
    175.768976227317  , 20.1147541286731  ,
    174.764584302758  , 20.672131655276   ,
    174.513490746954  , 21.245902298642   ,
    199.874435520754  , 21.9508201874443  ,
    199.372230707803  , 21.0655745472506  ,
    194.601364640814  , 20.4426233979296  ,
    199.62334196495   , 20.4918039038846  ,
    189.07720020507   , 19.5409834853072  ,
    199.62334196495   , 19.4590167458259  ,
    209.918372467684  , 17.5081965583815  ,
    224.231006071339  , 19.0163933478962  ,
    234.526036574073  , 19.475409862589   ,
    239.548013898209  , 18.2295087196126  ,
    253.609536244717  , 18.508197482914   ,
    258.631478166167  , 18.6065573391585  ,
    261.644671641186  , 19.0819669706144  ,
    250.596342769698  , 19.2786889944345  ,
    209.667296613224  , 20.1967220238201  ,
    207.909588618567  , 20.3770497752114  ,
    211.927191719488  , 20.6065580325578  ,
    224.231006071339  , 22.1639341723619  ,
    224.482081925799  , 21.7540981636242  ,
    230.006246361544  , 21.5081967895146  ,
    234.526036574073  , 20.9836066521037  ,
    238.0414171607    , 20.5573775266028  ,
    240.552388121425  , 20.6557373828472  ,
    250.09415565809   , 21.0655745472506  ,
    250.596342769698  , 21.5901646846616  ,
    250.345231512551  , 22.0491800436887  ,
    250.345231512551  , 22.2459020675089  ,
    255.11609757954   , 23.7704919737867  ,
    275.455118133234  , 21.4590162835596  ,
    274.952895618939  , 21.0655745472506  ,
    269.428731183195  , 20.9836066521037  ,
    268.173245702833  , 20.6885247720391  ,
    255.11609757954   , 20.4918039038846  ,
    269.679842440342  , 19.5081972517809  ,
    305.335835418219  , 19.1803279825244  ,
    304.833648306612  , 19.5573777577359  ,
    313.622082071836  , 20.2950818800645  ,
    300.313858094083  , 20.9344273018143  ,
    300.313858094083  , 21.8032798252448  ,
    300.313858094083  , 22.0327869269256  ,
    300.313858094083  , 22.6721311930097  ,
    344.507208982723  , 22.1475422112644  ,
    333.709991368382  , 22.7213116989648  ,
    365.850603759632  , 23.4918032104853  ,
    389.956045351726  , 22.1639341723619  ,
    384.93406802759   , 24.54098406314    ,
    398.242292005343  , 25.4918039038846  ,
    451.475223319041  , 25.0000005778328  ,
) |>
    mutate(
        workload = round(workload, 1),
        GE = round(GE, 3)
    )

model = lm(GE ~ log(workload), data = GE_data)

newdata <- data.frame(
    workload = floor(min(GE_data$workload)):ceiling(max(GE_data$workload))
)

pred <- predict(model, newdata = newdata, interval = "prediction", level = 0.90)

model_data <- cbind(newdata, pred) |>
    tibble() |>
    rename(GE = fit, conf.low = lwr, conf.high = upr)

ggplot(GE_data, aes(x = workload, y = GE)) +
    labs(
        title = "Gross Efficiency (GE) increases **logarithmically** with power output",
        subtitle = "Reproduced from Ettema & Lorås, 2009. Fig. 2d.",
        caption = str_glue(
            "**Data**: DOI: 10.1007/s00421-009-1008-7 | **Visuals**: {social_caption}"
        )
    ) +
    theme(
        plot.subtitle = ggtext::element_textbox(margin = margin(t = 4, b = 4))
    ) +
    scale_x_continuous(
        name = "Power Output (W)",
        expand = expansion(mult = 0.01)
    ) +
    scale_y_continuous(
        name = "Gross Efficiency (%)",
        expand = expansion(mult = 0.01)
    ) +
    geom_point(fill = "white", size = 2.4, shape = 21, stroke = 0.5) +
    geom_point(colour = "white", size = 1.6) +
    geom_ribbon(
        data = model_data,
        aes(ymin = conf.low, ymax = conf.high, colour = "GE", fill = "GE"),
        alpha = 0.2
    ) +
    geom_line(data = model_data, aes(colour = "GE")) +
    geom_errorbar(
        data = filter(model_data, workload == 250),
        aes(ymin = conf.low, ymax = conf.high),
        linewidth = 0.8, width = 20
    ) +
    annotate(
        "curve",
        x = 310, y = 15, xend = 250, yend = 18,
        curvature = -0.5, arrow = arrow(length = unit(5, "mm"))
    ) +
    annotate(
        "label", x = 350, y = 15, 
        label = "90% prediction intervals\n± 2.5% GE"
    ) +
    geom_point(data = filter(model_data, workload == 250), size = 3)